Stationary scattering from a nonlinear network 
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Transmission through a complex network of nonlinear one-dimensional leads is discussed by ex- 
tending the stationary scattering theory on quantum graphs to the nonlinear regime. We show that 
the existence of cycles inside the graph leads to a large number of sharp resonances that domi- 
nate scattering. The latter resonances are then shown to be extremely sensitive to the nonlinearity 
and display multi-stability and hysteresis. This work provides a framework for the study of light 
propagation in complex optical networks. 
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The study of quantum graphs has gained its popularity 
in recent years [l| not only because graphs emulate suc- 
cessfully complex mesoscopic and optical networks, but 
also because they manage to reproduce universal prop- 
erties (such as level statistics, transmission fluctuations 
and others) observed in generic quantum chaotic systems. 
Here we generalize quantum graph theory to the nonlin- 
ear domain. The theory will be applied in particular to 
show the effect of nonlinearity on transmission through 
networks of nonlinear fibers. Our model may also be 
used as a simple yet non-trivial model where the univer- 
sal properties derived from detailed numerical computa- 
tions of Bose-Einstein condensates in non-regular traps 
0-0| could be further investigated. 

Scattering is studied as a stationary process. The 
main finding is that the sharp resonances which dom- 
inate scattering in networks with complex connectivity 
lead to a dramatic amplification of the nonlinear effects: 
while non-resonant scattering hardly deviates from the 
predictions of the linear theory, tuning the parameters to 
a nearby resonance (without changing the incoming field 
intensity) brings the system into the nonlinear regime 
which is signaled by multi-stability and hysteresis. For 
this reason we revisit the theory of scattering in the lin- 
ear regime and demonstrate that sharp resonances with 
large amplification of the incoming wave inside the sys- 
tem are very frequent for graphs compared to other com- 
plex (chaotic) scattering systems. The origin of this effect 
can be related to the topology of the graph (existence of 
cycles) and leads to a power-law distribution for the am- 
plification. 



I. THE NONLINEAR SCHRODINGER 
EQUATION ON GRAPHS 

Consider a general metric graph which consists of V 
vertices connected by B internal bonds and N leads to 
infinity as illustrated in Fig. [TJ The bonds and leads 
will be collectively referred to as edges. The bonds are 
of finite lengths Lb, b = 1, • • • , B and are endowed with 
coordinates Xb e [0, Lb] (with a definite choice for the 
direction in which Xb increases). The semi-infinite leads 
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FIG. 1. (color online) A graph with V=4 vertices, B—6 bonds 
and N=2 leads. The incoming, reflected and transmitted 
waves are shown on the respective leads. 



have the coordinate xi € [0, oo) and xi — is at the ver- 
tex where the lead is attached. The wave function on 
the graph is a bounded piecewise continuous and differ- 
entiable function on the edges written collectively as 



B+N 
e=l 



(1) 



where ip e (x e ) is the wave functions on the edge e. While 
the model can describe far more general settings we re- 
strict ourselves in this exploratory work to the discus- 
sion of stationary scattering. The wave function on edge 
e satisfies the stationary nonlinear Schrodinger equation 
(NLSE) 



d 2 ^ e 
dx 2 



g e \lp e \ 2 ll) e = El() e 



(2) 



Here, g e is real nonlinear coupling parameter which we 
assume constant on each edge (but it may take different 
values on different edges). E is taken positive throughout 
this manuscript, E — k 2 and k reduces to the wave num- 
ber (propagation constant in fiber optics) in the linear 
case. Setting g e = on all edges will reduce our model 
to a standard (linear) quantum graph. Note that for ap- 
plications in fiber optics, g\^\ 2 >C E and the nonlinear 
term is a small perturbation of an otherwise linear wave 
equation. Then the stationary equation above describes 
the spatial evolution of the amplitude of a continuous 
wave beam rather than a wave envelope (for which one 
would have a non-stationary nonlinear Schrodinger equa- 
tion 
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A. The wave function on a single bond 

The solutions of © on a single bond can be obtained 
@ by writing i(>(x) = r(x)e l6 ^ (omitting the index e 
for the moment). Then, © is equivalent to two coupled 
ordinary differential equations 



dx 



and 



dx 



(3a) 



task. Since the purpose here is to display general fea- 
tures of wave propagation through a nonlinear network, 
we chose a "minimal" set of local matching conditions 
commonly used in the linear case: we require continuity 
and that the sum over all outgoing derivatives of wave 
functions on edges adjacent to j to be proportional to 
the common value 6a of the wave-function at the vertex, 
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These equations formally describe a classical particle in 
a central potential V(r) = (E/2)r 2 - (g/i)r 4 on the 2D 
plane where x now takes the role of time. The Hamil- 
tonian energy % and the angular momentum C are con- 
stants of motion. Note that the angular momentum pep 
reduces to the flux 



C — Im ip*dip/dx 



(4) 



carried by the wave function. For given values for % and 
C the full solution is obtained in the form of the integrals 
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which can be reduced to elliptic integrals @. 



B. Matching conditions at the vertices 

Two physical requirements guide our choice of the 
matching conditions at the vertices: continuity and the 
conservation of flux. Consider a vertex j with Vj ad- 
jacent edges and set x e — as the vertex coordinate 
on all the edges e emanating from j. Using the clas- 
sical point-particle analogue, continuity implies that at 
"time" x e = 0, all radii r e and angles 9 e assume the 
same values. Flux conservation is equivalent to angular 
momentum conservation, 



E 



dx e 



(7b) 



x c =0 



The constants Xj are arbitrary real parameters. In the 
classical particle picture, the imaginary part of this con- 
dition ensures conservation of angular momenta ([5]) and 
the real part can be expressed via 



I>(°) 



(8) 



e=l 



where p e is the radial momentum associated with the 
particle on the edge e, and rj is the radial coordinate 
at the vertex. When Vj = 2 the matching condition is 
equivalent to replacing the vertex by a S potential with 
strength Xj. 

In linear quantum graphs theory, it was found useful 
to express the matching conditions in terms of a vertex 
scattering matrix which connects the coefficients of 
incoming and outgoing waves Though in the nonlin- 
ear settings the lack of a superposition principle prohibits 
a decomposition into incoming and outgoing waves, the 
concept can be taken over formally. Defining 
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and collecting them in vectors a 11 

^a 1 1 n ^° ut '"', . . . , a 1 "^° ut '"'^ the matching conditions 
become 
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There is a large family of mathematically acceptable 
matching conditions which satisfy the latter requirements 
- e.g. all matching conditions that define a self-adjoint 
linear Schrodinger operator on the graph (see [l(|) satisfy 
flux conservation. The matching conditions appropriate 
for any particular experimental setting should in prin- 
ciple be derived ab initio, which is clearly a non-trivial 
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flux conservation 
of the vertex 
Im V*(0)# B (0)/da;e 



follows from the 
scattering matrix 



that J2e=i £e — becomes J2j l a 
- implying flux conservation. 
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assume that the vertex matching conditions (|7a|) and 
([7b|) are satisfied on all vertices. 

To finish the discussion of the vertex matching con- 
ditions we note that matching conditions for the time 
dependent NLSE on star graphs where discussed previ- 
ously i n f ill . Il2j | . The more relevant to the present work 
is Ref.JllJ where the authors treated rigorously the case 
t>j = 2 (see also 0)- For weak nonlinearity the vertex 
scattering matrix ()10bl) follows from their derivation. 



C. Scattering from a nonlinear network 

In the linear case the transport through a quantum 
graph with N leads can be described by an N X N uni- 
tary scattering matrix S(k) which connects incoming and 
outgoing amplitudes on the leads 



a 



out, leads 



S(k)z 



in. leads 
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The scattering matrix S(k) can be expressed explicitly 
in terms of the vertex scattering matrices <7«' , the bond 
lengths and the wave number k [13| in the form 



S(k) = p + T Q 
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Here T(k) is a diagonal 2B x 2B matrix with diagonal 
entries e lkLh that give the phase difference of a plane wave 
at the two ends of the bond b (each bond appears twice 
due to the two possible directions of a plane wave). The 
matrices p, cr lnt , r in , and r out are built up from the vertex 
scattering matrices . I.e. p is an N X N matrix that 
contains all direct scattering amplitudes (if all leads are 
attached to different vertices this is a diagonal matrix); 
the 2B x 2B matrix a- ln t contains all scattering amplitudes 
from one (directed) bond to another inside the graph; 
eventually r ln and r out are 2B x N and N x 2B matrices 
that contain scattering amplitudes from the leads into 
the bonds, and from the bonds out to the leads. They 
can be combined to one unitary matrix 



E = 



P Tout 
Tin (Tint 



(13) 



Unitarity of £ implies the unitarity of the scattering ma- 
trix S(k) The unitarity of S(k) ensures global flux con- 
servation 
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Moreover, in a linear system the scattering matrix S(k) 
is independent of the incoming amplitudes a m,loads . 

In the nonlinear case transport is described by a N- 
component nonlinear scattering function 
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(15) 



Flux conservation on each vertex implies that the scat- 
tering function conserves the norm 



\s(k, a 
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Though the general solution of the NLSE on each edge 
and the matching conditions are all known explicitly, 
it is generally not possible to solve the correspond- 
ing set of equations and obtain the scattering function 
s(k, a m ' loads ) in closed form. Therefore, we shall continue 
the discussion in the following section by presenting a nu- 
merical solution of a relevant example. 



II. RESONANT SCATTERING FROM A 
NONLINEAR NETWORK 

A. A simple examplary model 

We study the graph shown in Fig. [TJ The six bond 
lengths were chosen by a random number generator in 
the interval < L& < 1 and rationally independent 
within the numerical accuracy. We have kept the same 
set of lengths for all numerics that is presented in this 
manuscript (see [lij for the actual values). While we do 
not show results for different (random) choices we have 
checked that they lead to qualitatively equivalent results 
(the statistical properties that we will mention are also 
quantitatively equivalent). 

Two linear leads (L for 'left' and R for 'right') are 
attached at the vertices 1 and 2 respectively, with ql,r — 
0. A stationary wave with E = k 2 and intensity /; n = 



,in|2 



incident from the left lead is partially transmitted 



to the right lead and partially reflected 

M*l) =lfc» + ^refl =4" (e' tkXL + r(k, af) e^) 



=< t(k, a£) e ikXR 



(17a) 
(17b) 



Global gauge invariance implies that the phase of a 1 ^ 
can be chosen arbitrarily, so we take a™ = \/^in > 0. 
The reflection and transmission coefficients r(k, a 1 ^) and 
t{k, a™) will be computed as functions of k and I- m . In the 
linear case (g e — for all bonds) the reflection and trans- 
mission amplitudes are the matrix elements of a 2 x 2 scat- 
tering matrix S(k). In the nonlinear case they form the 
two components of the scattering (vector valued) func- 
tion s = (r(k),t(k)) T . Flux conservation implies 



\r(k,a}£)\ 2 + \t(k,a}£)\ 2 = l 



(18) 



In the present setting the importance of the nonlinear 
effects is controlled by L ln . The linear theory is obtained 
in the limit /;„ — > 0. However, the strength of the nonlin- 
earity will not be uniform as a function of k because the 
intensity inside the graph structure may vary strongly, 
especially near resonances as we will show below. 

Using the formalism described above, the solution of 
the scattering problem reduces to a finite set of non- 
linear equations in a high-dimensional space which re- 
quires rather substantial computer resources. The nu- 
merical complexity can be alleviated further by consid- 
ering the special case with just one nonlinear bond b, so 
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that gb = ±5 b j J . In the present simulation we choose 

b = (1,2) but we have also checked that other choices 
yield quantitatively similar behavior. 

B. Linear scattering: resonances and amplification 

The key to the understanding of the amplification of 
the nonlinear effects resides with the scattering in the 
linear regime. The upper panel of Fig. [5] shows the (lin- 
ear) reflection probability |r(fc,0)| 2 as a function of the 
wave number k for an interval of moderate wave num- 
bers where the typical length of a bond is about 14- 
16 wave lengths. Several resonances are clearly visible. 
Quantum graphs with incommensurate bond lengths are 
a paradigm of quantum chaotic scattering |13f where 
the statistics of resonances (i.e. of their location and 
their widths) is very close to the universal predictions 
of random-matrix theory. This implies that the width 
of the resonances is distributed over a broad range of 
values as illustrated by Fig. [5J The width of a single 
resonance is inversely proportional to the decay time of 
the corresponding resonant state. Narrow resonances are 
associated with waves which are trapped in the structure 
for long time, which is expressed in the stationary formal- 
ism by relatively large values of the wave function on the 
bonds. The lower panel in Fig. [5] shows the amplification 
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FIG. 2. (color online) Upper panel: reflection probability 
jr(fc,0)| 2 for the graph depicted in figure[T]in the linear limit. 
The arrow marks the narrow resonance which is discussed in 
the text. Lower panel: intensity amplification factor a(k) 
in logarithmic scale in the linear limit a™ —¥ on the bond 
connecting vertices 1 and 2 in figure [T] 

factor 

a(k) = ^ l^^)' 2 ^ (19) 

Lbhn 



for the bond (1,2) in the linear limit. While the in- 
tensity on the bond fluctuates around the incoming in- 
tensity, /in, there are also large peaks at narrow res- 
onances. For example, near the marked sharp reso- 
nance in Fig[2] the intensity on the bond is two orders 
of magnitude (« 320 times) higher than the intensity 
of the incoming beam. Over a larger spectral interval 
(0 < k < 20000) we found several other resonances 
with amplification factors a > 10 5 and a distribution 
P(a) = K^ 1 S(a — a(k))dk with a power law decay, 
P(a) ~ a~ s with s rj 2.85. The algebraic decay of this 
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FIG. 3. (color online) Double logarithmic plot of the numer- 
ically obtained distribution of the amplification factor on the 
nonlinear bond for the graph depicted in figure [1] 

distribution is a special feature of networks which would 
usually not be expected in other chaotic scattering mod- 
els such as scattering through a chaotic quantum dot. To 
explain the difference between a network and a more gen- 
eral chaotic scattering system let us go back to Eq. p^|) 
which describes scattering through the network. Simi- 
lar equations have been used to model chaotic scattering 
(e.g. by an average over the matrix E from Eq. Q13p). 
In either case the factor (1 — T(k)cr- m t) in the second 
term is responsible for the resonances and one may ex- 
pect large amplification factors whenever T(k)ai n t has 
an eigenvalue near unity (an exact eigenvalue unity in- 
dicates the existence of a bound state which is confined 
to the bonds). In any generic model of chaotic scatter- 
ing unimodular eigenvalues of the subunitary matrix <7i n t 
(which acts on a vector of 2B coefficients, one for each 
directed bond) are strongly suppressed. For other known 
examples of resonant scattering in one-dimensional non- 
linear Schrodinger systems [5j-l7| , the equivalent of (Tint is 
one number with modulus smaller than 1. In all these 
models high amplification factors are either cut-off or ex- 
tremely rare. However for networks with the standard 
vertex matching conditions (|7aj) and (|Tb[) (and Xj = 
see below for Xj ^ 0) the situation is drastically differ- 
ent as every cycle created from the bonds of the network 
supports an eigenvalue unity of dint (moreover cycles of 
even length support eigenvalues minus unity). In fact, 
let us consider a cycle that consists of three bonds b\, 
62, and 63, then corresponding eigenvector a of <7i n t with 
unit eigenvalue vanishes on all directed bonds that do not 



5 



belong to the cycle, and has values ±1 on the directed 
bonds that belong to the cycle (the two signs correspond 
to two different directions to go through the cycle). For 
rationally dependent bond lengths this implies that one 
may chose k such that 

e ikL bl _ e ikL b2 _ £ ikL b3 _ ^ ^0) 

which shows the existence of embedded bound states in 
the continuum of scattering states. The construction is 
equivalent to that of perfectly scarred states (see (l6| ) 
- these scarred states vanish exactly on the vertex and 
for each bond on the cycle the bond length is an integer 
multiple of the wave length 2it/k. For incommensurate 
bond lengths, there are no perfectly scarred states on 
the cycle as the condition (|2T)|) can never be met exactly. 
However the mapping k — > (e lfcif, i , e b i , e lkLb z ) is an er- 
godic flow on a 3-dimensional torus, one thus finds values 
for the wave number k which approximate condition (|20[) 
to arbitrary precision. In the exemplary model we have 
used for our calculation the graph structure contains two 
independent cycles that both contain the nonlinear bond. 
The two cycles can be chosen such that they consist of 
three bonds. 

In the above discussion we have assumed the the ver- 
tex potentials Xj vanish. The vertex scattering matrices 
(llObp show however that for sufficiently high wave num- 
ber k these potentials are not relevant. Note also, that 
vertex matching conditions which are entirely different 
from (|7ap and (|7b[) do not necessarily have a similar dis- 
tribution of narrow resonances. 



C. Implications for nonlinear scattering: 
multistability 

The strength of nonlinearity on the nonlinear bond 
may be measured by the effective parameter 

v = M]^^Mj^ lMxe) f dXe . (21) 

For a fixed incoming intensity /;„, and with \g e \ = 1 the 
effective nonlinearity will proportional to the amplifica- 
tion factor 

u(E) = ^a(E)l m . (22) 

Even if the incoming intensity is too low to induce notice- 
able nonlinear effects off resonance, at narrow resonances 
the high fields on the bonds are expected to behave in 
a nonlinear way. This qualitative picture is supported 
by the numerical simulations. For incoming intensities 
up to /in ~ 0.005 the nonlinearity is either not relevant 
at all or can be taken into account as a perturbation for 
almost the entire fc-spectrum. However, near the marked 
resonance the amplification by two orders of magnitude is 
sufficient to give rise to strong nonlinear effects that can- 
not be described perturbatively. This is shown in figure 
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FIG. 4. (color online) Reflection probability near the narrow 
resonance marked in FigO The central (blue) curve is the res- 
onance in the linear limit. The (green/red) curves left/right 
from the central one correspond to 8 increasing values of the 
incoming intensity \a^\ 2 (in equal steps from 0.0006 to 0.0048) 
in the attractive/repulsive case (g — 1/ g — —1). The inset 
shows the corresponding amplification factor. 

Fig. U which resolves the narrow resonance in the reflec- 
tion probability (and the amplification factor) for vari- 
ous values of the incoming intensity. In the attractive 
(repulsive) case the resonance moves to the left (right) 
as I- m increases. However some parts of the curve move 
faster than others which eventually leads to a multivalued 
dependency above the critical value / C rit- This implies 
multistability - an experiment would show hysteresis. 
In both the attractive and the repulsive cases the crit- 
ical incoming intensity where multistability set in is near 
J crit w 0.002. For this incoming intensity the strength of 
nonlinearity inside the graph is typically (i.e. away from 
the resonance) on the order of v tTP ~ 5 x 10~ 8 - at the 
resonance it is however ^ rcs = 1.5 x 10~ 5 . These findings 
are similar to previous work on nonlinear resonant scat- 
tering from quantum dots and from one-dimensional 
structures [5|-0] . Our model generalises the latter results 
by allowing for additional topological complexity. 



D. Application in nonlinear fiber optics 

The numerically found power-law distribution of the 
amplification factor a can be expected to be a generic 
feature (at least for similar types of vertex matching 
conditions). This implies that by tuning the parameters 
of an experiment to a sufficiently narrow resonance one 
may find arbitrarily high amplification factors and thus 
the nonlinearity effects of multistability and hysteresis 
may be observed at considerably lower incoming veloci- 
ties then in our example. Nonetheless let us translate the 
important parameters of our model to a fiber- network ex- 
perimental setup. 

For a CW optical beam propagating in a single mode 
telecom fiber with the linear refraction index no = 1.5, 
nonlinear Kerr coefficient ri2 — 2.4 x 10~ 16 cm 2 /W, and 
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effective mode area S e g = 50/im 2 operating at the tele- 
com wavelength A = 2-kc/uo = 1.55/i the strength of the 
nonlinearity (on the nonlinear bond) can be estimated as 



8n w 2 
c 2 f3 2 S cS 



(23) 



Here, (3 is the propagation constant and P avc is the bond 
averaged power \ip\ 2 of the beam in watts. For the pur- 
poses of the current estimate one can neglect mode dis- 
persion of the fiber and assume j3 ~ no(2-7r/A). Then 
from (f2"3")l it follows that in order to achieve the resonance 
level of nonlinearity v Tes = 1.5xl0~ 5 from the numerical 
example considered above at the telecom wavelength A, 
the required power level must be as high as P avc ~ 6kW 
which is above the thresholds for stimulated Raman and 
Brillouin scattering for the fiber length of a few meters. 
This can be offset in several ways: one can consider mode 
dispersion and operate at lower frequencies so that the 
effective mode index /3c/cj is lower; or, reduce the effec- 
tive mode area by a factor up to 10 by changing the size 
of the fiber core or else one could use highly nonlinear 
non-silica fibers where the value of the nonlinear coeffi- 
cient n.2 can be enhanced up to two orders of magnitude 
Thus, using the estimates based on the present sim- 
ulation, one can expect that nonlinear effects will appear 
at power levels of 1 -j- 10W. 

In an experiment varying the wave number in a con- 
trolled way may not be feasible - so let us mention that 



one may equivalently change the lengths of the edges in a 
controlled way by slowly varying the temperature. Let us 
also mention that multistability in scattering from non- 
linear cristals has been observed experimentally 



III. CONCLUSION 

To conclude, the theory presented here shows how 
the interplay between complex topology and nonlinear- 
ity gives rise to a pronounced amplification of nonlinear 
effects. We would like to stress that while our model 
is highly idealized, the strong amplification of intensity 
near narrow resonances is a universal effect that can be 
expected in any linear complex network. Any co-existing 
nonlinearity that may be negligible off resonance will be 
drastically amplified at a resonance. We believe that the 
latter effect can be observed in actual experiments with 
interconnected optical fibers even though the model itself 
may need further adjustment to fit the details of such an 
experiment. Moreover we believe that NLSE on metric 
graphs as presented here will be a very useful paradigm 
system where the interplay between topology and non- 
linearity can be studied qualitatively. 
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